Load identification method for reciprocating machinery based on information entropy and envelope features of axis trajectory of piston rod

ABSTRACT

A load identification method for reciprocating machinery based on information entropy and envelope features of an axis trajectory of a piston rod. According to the present disclosure, firstly, the position of an axial center is calculated according to a triangle similarity theorem to obtain an axial center distribution; secondly, features are extracted from the axial center distribution of the piston rod by means of an improved envelope method for discrete points as well as an information entropy evaluation method; thirdly, a dimensionality reduction is carried out on the features by means of manifold learning to form a set of sensitive features of the load; and finally, a neural network is trained to obtain a load identification classifier to fulfill automatic identification on the operating load of the reciprocating machinery. The advantages of the present disclosure are verified by means of actual data of a piston rod of a reciprocating compressor.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority to Chinese Application Serial No.201911079879.2, filed on Nov. 7, 2019, the content of which isincorporated herein by reference in its entirety.

TECHNICAL FIELD

The present disclosure relates to a load identification method forreciprocating machinery.

BACKGROUND

A load variation typically refers to changes of dynamic characteristicsof mechanical structures, which lead to an influence on fault featuresof vibration signals, displacement signals, and the like. It has alwaysbeen difficult to implement fault monitoring and diagnosis in variableload conditions. Piston rods as movable key parts of reciprocatingmachinery are prone to causing loosening, cracking, or even breaking tofasteners. There have been many research reports on the fault monitoringand diagnosis of the reciprocating machinery as well as research reportson axis trajectories of the piston rods. For example, an acousticemission technology is adopted to perform on-line monitoring on thepiston rods to pre-warn accidents. A method for fault diagnosis andanalysis based on the axis trajectories of the piston rods in the Xdirection and Y direction can fulfill an early warning of potentialfaults on the piston rods and piston assemblies of the reciprocatingmachinery. A harmonic wavelet is used to extract vibration energy,natural frequencies, areas of trajectory envelopes, and other featuresbased on the axis trajectories of the piston rods for fault diagnosis.

At present, there are few reports on the extraction of motion featuresof the piston rods in the variable load conditions. When a change to theload conditions and faults occur simultaneously, it is necessary todetermine why the load conditions of the piston rods change, and thefeatures of influences on loads and the faults should be deeplyresearched to prevent the faults from being falsely determined.Accordingly, it can be valuable to research transient motion features ofthe piston rods in the variable load conditions. In view of this, thepresent disclosure provides a method for extracting information entropyand envelope features of a discrete point distribution contour based onan axis trajectory of a piston rod, which extracts the features of theaxis trajectory of the piston rod in different load conditions andestablishes a set of sensitive feature parameters of a load as well as aload identification model by means of training.

SUMMARY

The objective of the present disclosure is to provide a simple andeffective load identification method for reciprocating machinery, whichextracts information entropy and envelope features based on data of anaxis trajectory to establish a set of sensitive features of the load forload identification on reciprocating machinery. The present disclosurehas simple calculation, high adaptability, high accuracy inidentification, and the like.

The objective of the present disclosure is fulfilled through thefollowing technical solution: firstly, the position of an axial centeris calculated based on settlement data and deflection data of a pistonrod; secondly, an envelope feature of an axial center distribution isextracted by means of an improved envelope method for a discrete pointdistribution contour, then an information entropy feature of the axialcenter distribution is calculated, and an initial feature set is formedby the envelope feature and the information entropy feature; andthirdly, sensitive features of the load are extracted from the initialfeature set by means of manifold learning to form a set of the sensitivefeatures of the load, and a neural network is trained by means of theset of the sensitive features of the load to obtain an identificationclassifier.

A load identification method for reciprocating machinery based oninformation entropy and envelope features of an axis trajectory of apiston rod includes the following steps:

step 1. setting different load conditions Load={0, d, 2d, 3d, . . . ,wd}, w=0, 1, 2, . . . , where d represents a load gradient, and thenumber of the load conditions is (w+1) in total; respectively acquiring,by an on-line monitoring system of reciprocating machinery, an originaldeflection displacement X_(m)={x₁, x₂, x₃, . . . , y_(m)} and originalsettlement displacement Y_(m)={y₁, y₂, y₃, . . . , y_(m)} of a pistonrod in a corresponding load condition through an eddy currentdisplacement sensor (deflection sensor) in a horizontal direction and aneddy current displacement sensor (settlement sensor) in a verticaldirection to obtain an original data set XY_(n)={(X_(m),Y_(m))₁ ^(T),(X_(m),Y_(m))₂ ^(T), . . . , (X_(m),Y_(m))_(n) ^(T)}^(T), where mrepresents the number of sampling points, and n represents the number ofdata groups;

step 2. removing average values of an original signal X_(m) and anoriginal signal Y_(m) by means of formula (1) to obtain X′_(m)={x′₁,x′₂, x′₃, . . . , x′_(m)} and Y′_(m)={y′₁, y′₂, y′₃, . . . , y′_(m)},where the original data set is turned to XY′_(n)={(X′_(m),Y′_(m))₁ ^(T),(X′_(m),X′_(m))₂ ^(T), . . . , (X′_(m),Y′_(m))_(n) ^(T)}^(T)

$\begin{matrix}{{{F_{m}^{\prime}(j)} = {{F_{m}(i)} - {\frac{1}{m}{\sum\limits_{i = 1}^{m}{F_{m}(i)}}}}}{i,{j = 1},2,\ldots\mspace{14mu},m}} & (1)\end{matrix}$

where, in formula (1), F_(m) represents the original deflection orsettlement displacement of the piston rod, and F′_(m) represents thedeflection or settlement displacement, obtained after the average valuesare removed, of the piston rod; and

setting a horizontal direction measured by a deflection sensor as anX-axis and a vertical direction measured by a settlement sensor as aY-axis to build a plane-rectangular coordinate system, where if aposition of an axial center of the piston rod at an initial time isdenoted by O₀(a₀, b₀), the position of the axial center of the pistonrod at another time is denoted by O_(m)(a_(m),b_(m)), and a radius ofthe piston rod is denoted by R, a point of intersection between acircumference of the piston rod and the X-axis at this time isJ_(X)(R+x′_(m),0), and a point of intersection between the circumferenceof the piston rod and the Y-axis at the time is J_(Y)(0,R+y′_(m));setting an included angle between the X-axis and a line connecting thepoint O_(m) to the point J_(X) as θ and an included angle between a lineconnecting the point O_(m) to the point J_(Y) and a straight line,parallel to the X-axis, on which the point O_(m) is located as φ toderive formula (2) and formula (3) according to a triangle similaritytheorem; and solving, by means of formula (2) in combination withformula (3), the position O_(m)(a_(m),b_(m)) of the axial center of thepiston rod at different times to form an axial center distribution setO={O₁(a₁,b₁), O₂(a₂,b₂), O₃(a₃,b₃), . . . , O_(m)(a_(m),b_(m))}

$\begin{matrix}\{ \begin{matrix}{\frac{R + {X_{m}^{\prime}(j)} + a_{m}}{X_{m}^{\prime}(j)} = \frac{R}{{{X_{m}^{\prime}(j)}/\cos}\;\theta}} \\{\theta = {\arctan\frac{b_{m}}{R + {X_{m}^{\prime}(j)} + a_{m}}}}\end{matrix}  & (2) \\\{ \begin{matrix}{\frac{R + {Y_{m}^{\prime}(j)} - b_{m}}{Y_{m}^{\prime}(j)} = \frac{R}{{{Y_{m}^{\prime}(j)}/\sin}\;\varphi}} \\{\varphi = {\arctan\frac{R + {Y_{m}^{\prime}(j)} - b_{m}}{a_{m}}}}\end{matrix}  & (3)\end{matrix}$

where, in formula (2) and formula (3), j=1, 2, 3, . . . , m

step 3. calculating an envelope feature B_(ao) of the axial centerdistribution O={O₁(a₁,b₁), O₂(a₂,b₂), O₃(a₃,b₃), . . . ,O_(m)(a_(m),b_(m))} by means of an improved envelope method for adiscrete point distribution contour, where the improved envelope methodfor the discrete point distribution contour particularly includes thefollowing steps:

step 3.1. determining, according to the axial center distribution O,four limit points by seeking a minimum point a_(l) and a maximum pointa_(r) in the horizontal direction X as well as a minimum point b_(d) anda maximum point b_(u) in the vertical direction Y, where the four limitpoints are respectively denoted by O_(l)(a_(l),b_(l)),O_(r)(a_(r),b_(r)), O_(d)(a_(d),b_(d)), O_(u)(a_(u),b_(u)), an inside ofa quadrangle formed by the four limit points is counted as an internalside, and an outside of the quadrangle is counted as an external side;

step 3.2. extracting a convex envelope of an axial center distributioncontour at a minimum slope with the foregoing limit points as startingpoints by anticlockwise traversing all over the positions of the axialcenter at all times; and

calculating a convex envelope between a limit point O_(d) and a limitpoint O_(r) through a method including the following steps;

(1) setting a line connecting the point O_(d) to the point O_(r) as L₁,where a slope α₁ of the line is expressed as:

$\begin{matrix}{\alpha_{1} = \frac{b_{d} - b_{r}}{a_{d} - a_{r}}} & (4)\end{matrix}$

(2) if a set of all common axial center points on an external side ofthe line L₁ is assumed as P={p₁(a₁,b₁), p₂(a₂,b₂), p₃(a₃,b₃), . . . },calculating a slope K={β₁, β₂, β₃, . . . } of a line connecting thepoint O_(d) to any point in P; if there are multiple points with a sameslope, calculating a distance D={dis₁, dis₂, dis₃, . . . } between thecorresponding point and the point O_(d); and seeking a pointp′(a_(p),b_(p)) in the set P, and enabling p′ to meet the followingformula:β′=min K,β′≤α ₁, and dis′=max D  (5)

where, in formula (5), β′ represents a slope of a line connecting thepoint p′ to the point O_(d), and dis′ represents a distance between thepoint p′ and the point O_(d); and

the point p′ is a convex envelope point;

(3) replacing the point O_(d) with the convex envelope point p′, settinga line connecting the point p′ to the point O_(r) as L′₁ and a slope ofthe line L′₁ as α′₁, and then carrying out a next iteration to seek anew convex envelope point;

(4) repeating step (1), step (2), and step (3); and when a distancebetween the new convex envelope point and the point O_(r) is 0, stoppingthe iteration to obtain a convex envelope B′_(dr)={p′₁, p′₂, p′₃, . . .} of an axial center distribution contour between the limit point O_(d)and the limit point O_(r); and

calculating a convex envelope B′_(ru) between the limit point O_(r) anda limit point O_(u), a convex envelope B′_(ul) between the limit pointO_(u) and a limit point O_(l), and a convex envelope B′_(ld) between thelimit point O_(l) and the limit point O_(d) to obtain a setB_(tu)={O_(d), B′_(dr), O_(r), B′_(ru), O_(u), B′_(ul), O_(l), B′_(ld)}of all convex envelope points in the axial center distribution;

(5) calculating an area S1 of the convex envelope of the axial centerdistribution contour by means of formula (6) below:

$\begin{matrix}{{S\; 1} = {\frac{1}{2}{\sum\limits_{e = 1}^{c\; 1}( {{a_{pe}*b_{p{({e + 1})}}} - {a_{p{({e + 1})}}*b_{pe}}} )}}} & (6)\end{matrix}$

where, in formula (6), c1 represents the number of the convex envelopepoints;

step 3.3. calculating a concave envelope of the axial centerdistribution according to the convex envelope obtained in step 3.2; and

calculating a concave envelope between the limit point O_(d) and thelimit point O_(r) by successively anticlockwise carrying out thefollowing calculation on two continuous convex envelope pointsp′₁(a_(p1),b_(p1)) and p′₂(a_(p2),b_(p2)) in the convex envelopeB′_(dr)={p′₁, p′₂, p′₃, . . . }:

(1) setting a line connecting the point p′₁(a_(p1),b_(p1)) to the pointp′₂(a_(p2),b_(p2)) as L₂, where a slope α₂ of the line L′ is expressedas:

$\begin{matrix}{\alpha_{2} = \frac{b_{p\; 2} - b_{p\; 1}}{a_{p\; 2} - a_{p\; 1}}} & (7)\end{matrix}$

(2) if a set of all common axial center points on an internal side ofthe line L₂ is assumed as Q={q₁(a₁,b₁),q₂(a₂,b₂),q₃(a₃,b₃), . . . },calculating a slope K′={β′₁, β′₂, β′₃, . . . } of a line connecting thepoint p′₁ to any point in Q; if there are multiple points with a sameslope, calculating a distance D′={dis′₁, dis′₂, dis′₃, . . . } betweenthe corresponding point and the point p′₁; and seeking a pointq′(a_(q),b_(q)) in the set Q, and enabling q′ to meet the followingformula:β″=min K′,β″≥α ₂, and dis″=max D′  (8)

where, in formula (8), β″ represents a slope of a line connecting thepoint q′ to the point p′₁, and dis″ represents a distance between thepoint q′ and the point p′₁; and

the point q′ is a concave envelope point;

(3) replacing the point p′₂ with the concave envelope point q′, settinga line connecting the point p′₁ to the point q′ as L′₂ and a slope ofthe line L′₂ as α₂, and then carrying out a next iteration to seek a newconcave envelope point;

(4) repeating step (1), step (2), and step (3); and when a distancebetween the new concave envelope point and the point p′₁ is not greaterthan M, stopping the iteration to obtain a concave envelopeB″_(dr)={q′₁, q′₂, q′₃, . . . } of the axial center distribution contourbetween the limit point O_(d) and the limit point O_(r); where,

an initial value of M indicates an average distance between two adjacentaxial center points, and M is calculated by means of formula (9) incombination with formula (10):

$\begin{matrix}{\mspace{20mu}{M = \frac{2M^{\prime}}{m( {m - 1} )}}} & (9) \\{M^{\prime} = {{\sum\limits_{i = 2}^{m}{\sqrt{( {a_{1} - a_{i}} )^{2} + ( {b_{1} - b_{i}} )^{2}}}} + {\sum\limits_{i = 3}^{m}{\sqrt{( {a_{2} - a_{i}} )^{2} + ( {b_{1} - b_{i}} )^{2}}}} + \ldots + {\sum\limits_{i = m}^{m}{\sqrt{( {a_{m - 1} - a_{i}} )^{2} + ( {b_{m - 1} - b_{i}} )^{2}}}}}} & (10)\end{matrix}$

calculating a concave envelope B″_(ru) between the limit point O_(r) anda limit point O_(u), a concave envelope B″_(ul) between the limit pointO_(u) and a limit point O_(l), and a concave envelope B″_(ld) betweenthe limit point O_(l) and the limit point O_(d) to obtain a setB_(ao)={O_(d), B″_(dr), O_(r), B″_(ru), O_(u), B″_(ul), B″_(ld)} of allconcave envelope points in the axial center distribution;

(5) calculating an area S2 of the concave envelope of the axial centerdistribution contour by means of formula (11) below:

$\begin{matrix}{{S2} = {\frac{1}{2}{\sum\limits_{g = 1}^{c2}( {{a_{qg}*b_{q{({g + 1})}}} - {a_{q{({g + 1})}}*b_{qg}}} )}}} & (11)\end{matrix}$

where, in formula (11), c2 represents the number of the concave envelopepoints;

step 3.4. determining whether or not the set B_(ao), obtained in step3.3, of the concave envelope points is the envelope feature of the axialcenter distribution of the piston rod;

(1) calculating a relative error E of S2 and S1 by means of formula (12)till E is less than or equal to 5%, where the set B_(ao), obtained instep 3.3, of the concave envelope points is the envelope feature of theaxial center distribution of the piston rod after the calculation isstopped;

$\begin{matrix}{{E = {\frac{{{S\; 2} - {S\; 1}}}{S1} \times 100\%}};} & (12)\end{matrix}$and

(2) in a case where the distance M is reduced by 50% when E is greaterthan 5%, replacing S1 with S2, repeating step 3.3 to obtain a new setB′_(ao) of the concave envelope points as well as an area S2′ of theconcave envelope; and repeatedly calculating a relative error E′ of S2′and S2 by means of formula (13) till E′ is less than or equal to 5%,where the iteration is stopped at this moment, and the set B′_(ao)obtained in Step 3.3 is the envelope feature of the axial centerdistribution of the piston rod during the last iteration;

$\begin{matrix}{E^{\prime} = {\frac{{{S\; 2^{\prime}} - {S\; 2}}}{S\; 2} \times 100\%}} & (13)\end{matrix}$

step 4. calculating an information entropy feature of the axial centerdistribution O: calculating, by means of formula (14), arithmetic squareroots of coordinates of all the points in the axial center distributionO to obtain a set S_(m)={s₁, s₂, s₃, . . . , s_(m)}, and thencalculating the information entropy feature Sh of the axial centerdistribution O by means of formula (15), where an initial feature setT={B_(ao),Sh} is formed by the information entropy feature Sh and theenvelope feature;

$\begin{matrix}{{{S_{m}( i^{\prime} )} = {\sqrt{a_{k}^{2} + b_{k}^{2}}}},k,{i^{\prime} = 1},2,\ldots\mspace{14mu},m} & (14) \\{{Sh} = {- {\sum\limits_{i^{\prime} = 1}^{m}( {\frac{{S_{m}( i^{\prime} )}}{\sum\limits_{i^{\prime} = 1}^{m}{{S_{m}( i^{\prime} )}}}*{\log_{2}( \frac{{S_{m}( i^{\prime} )}}{\sum\limits_{i^{\prime} = 1}^{m}{{S_{m}( i^{\prime} )}}} )}} )}}} & (15)\end{matrix}$

step 5. carrying out, by means of t-distributed stochastic neighborembedding (T-SNE), an unsupervised dimensionality reduction on theinitial feature set to extract sensitive features of a load; andassuming that the initial feature set T includes 1*Col-dimensionalfeatures, given perplexity is 30, and a given learning rate is 1e-5,setting a label as Labels={0, 1, 2, . . . , w}, and then inputting theinitial feature set T to a T-SNE algorithm for the unsuperviseddimensionality reduction in the (w+1) load conditions to obtain a setT′={t₁,t₂} of the sensitive features of a 1*2-dimensional load; and

step 6. firstly, sorting data acquired by the on-line monitoring systemin the (w+1) load conditions to form a training set and a test set;secondly, processing the data in the training set and the test setthrough the above steps to obtain a final training set Train_T′ and afinal test set Test_T′; and thirdly, setting, according to differentreciprocating machinery, the number of neurons of a back-propagation(BP) neural network as 20-30, a learning rate as 0.0005-0.001, trainingaccuracy as 0.0001-0.0005, and maximum iterations as 70-100, theninputting the data set Train_T′ to the BP neural network for training toobtain a classifier capable of distinguishing the (w+1) load conditionsof the reciprocating machinery, and testing the classifier of the BPneural network by means of the test set Test_T′.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a flow chart of a method of the present disclosure;

FIG. 2 is a schematic diagram showing the position of an axial center;

FIG. 3 shows a settlement waveform and deflection waveform of a pistonrod of a reciprocating compressor;

FIG. 4 shows an axial center distribution;

FIG. 5 shows an envelope feature calculated by means of an improvedmethod;

FIG. 6 shows sensitive features of a load; and

FIG. 7 shows an envelope feature calculated by means of a traditionalmethod.

DETAILED DESCRIPTION

For the sake of better understanding of the technical solution of thepresent disclosure, the specific embodiments of the present disclosureare further expounded below with reference to data of a piston rod of areciprocating compressor and the accompanying drawings.

Step 1. Acquire, by an on-line monitoring system of a reciprocatingcompressor, deflection data X_(m) and settlement data Y_(m) of a pistonrod in load conditions of 0%, 20%, 50%, 60%, 70%, 80%, 90%, and 100% toform an original data set XY_(n)={(X_(m),Y_(m))₁ ^(T), (X_(m),Y_(m))₂^(T), . . . , (X_(m),Y_(m))_(n) ^(T)}^(T), where n represents the numberof data groups and is equal to 500, waveforms are shown in FIG. 3,original data of eight loads fall into a training set and a test set,and 400 data groups in the training set and 100 data groups in the testset are set for each load;

Step 2. Remove average values of an original signal X_(m) and anoriginal signal Y_(m) to obtain X′_(m)={x′₁, x′₂, x′₃, . . . , x′_(m)}and Y′_(m)={y′₁, y′₂, y′₃, . . . , y′_(m)}, and solve, according to atriangle similarity theorem, the position of an axial center of thepiston rod to obtain an axial center distribution set O={O₁(a₁,b₁),O₂(a₂,b₂), O₃(a₃,b₃), . . . , O_(m)(a_(m),b_(m))}, where an axial centerdistribution is shown in FIG. 4;

Step 3. Calculate an envelope feature B_(ao) of the axial centerdistribution O by means of an improved envelope method for a discretepoint distribution contour, as shown in FIG. 5;

Step 4. Calculate an information entropy feature Sh of the axial centerdistribution, where an initial feature set T={B_(ao),Sh} is formed bythe information entropy feature and the envelope feature;

Step 5. In a case where given perplexity is 30 and a given learning rateis 1e-5, extract sensitive features of the loads from T by means ofmanifold learning to form a set T′={t₁,t₂} of the sensitive features ofthe loads, as shown in FIG. 6; and

Step 6. Process the data from the training set and the test set throughStep 2 to Step 5 to respectively obtain a final training feature setTrain_T′ and a final test feature set Test_T′; and set the number ofneurons of a back-propagation (BP) neural network as 30, a learning rateas 0.001, training accuracy as 0.0001, and maximum iterations as 100,input the data set Train_T′ to the BP neural network for training toobtain a classifier capable of distinguishing 8 load conditions of thereciprocating compressor, test the classifier of the BP neural networkby means of the test set Test_T′, and compare the improved envelopemethod with a traditional envelope method for extracting the envelopefeature (as shown in FIG. 7); where, a result is shown in Table 1.

TABLE 1 identification accuracy of the neural network (100 sets of testdata/load conditions) Overall identification Load/% 0 20 50 60 70 80 90100 accuracy/% Traditional 93 100 77 77 98 96 63 51 81.88 envelopemethod + information entropy Improved envelope 99 100 97 100 100 100 4774 89.63 method + information entropy

What is claimed is:
 1. A load identification method for reciprocatingmachinery based on information entropy and envelope features of an axistrajectory of a piston rod, comprising the following steps: step 1.setting different load conditions Load={0, d, 2d, 3d, . . . , wd}, w=0,1, 2, . . . , wherein d represents a load gradient, and the number ofthe load conditions is (w+1) in total; respectively acquiring, by anon-line monitoring system of reciprocating machinery, an originaldeflection displacement X_(m)={x₁, x₂, x₃, . . . , x_(m)} and originalsettlement displacement Y_(m) {Y₁, Y₂, Y₃, . . . , y_(m)} of a pistonrod in a corresponding load condition through an eddy currentdisplacement sensor in a horizontal direction and an eddy currentdisplacement sensor in a vertical direction to obtain an original dataset XY_(n)={(X_(m),Y_(m))₁ ^(T) (X_(m),Y_(m))₂ ^(T), . . . ,(X_(m),Y_(m))_(n) ^(T)}^(T), wherein m represents the number of samplingpoints, and n represents the number of data groups; step
 2. removingaverage values of an original signal X_(m) in and an original signalY_(m) by means of formula (1) to obtain X′_(m)={x′₁, x′₂, x′₃, . . . ,x′_(m)} and Y′_(m)={y′₁, y′₂, y′₃, . . . , y′_(m)}, wherein the originaldata set is turned to XY′_(n)={(X′_(m),Y′_(m))₁ ^(T), (X′_(m),Y′_(m))₂^(T), . . . , (X′_(m),Y′_(m))_(n) ^(T)}^(T) $\begin{matrix}{{{F_{m}^{\prime}(j)} = {{F_{m}(i)} - {\frac{1}{m}{\sum\limits_{i = 1}^{m}{F_{m}(i)}}}}}{i,{j = 1},2,\ldots\mspace{14mu},m}} & (1)\end{matrix}$ wherein, in formula (1), F_(m) represents the originaldeflection or settlement displacement of the piston rod, and F′_(m)represents the deflection or settlement displacement, obtained after theaverage values are removed, of the piston rod; and setting a horizontaldirection measured by a deflection sensor as an X-axis and a verticaldirection measured by a settlement sensor as a Y-axis to build aplane-rectangular coordinate system, wherein if a position of an axialcenter of the piston rod at an initial time is denoted by O₀(a₀,b₀), theposition of the axial center of the piston rod at another time isdenoted by O_(m)(a_(m),b_(m)), and a radius of the piston rod is denotedby R, a point of intersection between a circumference of the piston rodand the X-axis at this time is J_(x) (R+x′_(m),0), and a point ofintersection between the circumference of the piston rod and the Y-axisat the time is J_(Y) (0,R+y′_(m)); setting an included angle between theX-axis and a line connecting the point O_(m) to the point J_(x) as θ andan included angle between a line connecting the point O_(m) to the pointJ_(y) and a straight line, parallel to the X-axis, on which the pointO_(m) is located as φ to derive formula (2) and formula (3) according toa triangle similarity theorem; and solving, by means of formula (2) incombination with formula (3), the position O_(m) (a_(m), b_(m)) of theaxial center of the piston rod at different times to form an axialcenter distribution set $\begin{matrix}{{O = \{ {{O_{1}( {a_{1},b_{1}} )},{O_{2}( {a_{2},b_{2}} )},{O_{3}( {a_{3},\ b_{3}} )},\ldots\mspace{14mu},{O_{m}( {a_{m},b_{m}} )}} \}}\{ \begin{matrix}{\frac{R + {X_{m}^{\prime}(j)} + a_{m}}{X_{m}^{\prime}(j)} = \frac{R}{{{X_{m}^{\prime}(j)}/\cos}\theta}} \\{\theta = {\arctan\frac{b_{m}}{R + {X_{m}^{\prime}(j)} + a_{m}}}}\end{matrix} } & (2) \\\{ \begin{matrix}{\frac{R + {Y_{m}^{\prime}(j)} - b_{m}}{Y_{m}^{\prime}(j)} = \frac{R}{{{Y_{m}^{\prime}(j)}/\sin}\;\varphi}} \\{\varphi = {\arctan\frac{R + {Y_{m}^{\prime}(j)} - b_{m}}{a_{m}}}}\end{matrix}  & (3)\end{matrix}$ wherein, in formula (2) and formula (3), j=1, 2, 3, . . ., m step
 3. calculating an envelope feature B_(ao) of the axial centerdistribution O={O₁ (a₁, b₁), O₂ (a₂, b₂), . . . , O₃ (a₃, b₃), . . . ,O_(m) (a_(m), b_(m))} by means of an improved envelope method for adiscrete point distribution contour, wherein the improved envelopemethod for the discrete point distribution contour particularlycomprises the following steps: step 3.1. determining, according to theaxial center distribution O, four limit points by seeking a minimumpoint a_(l) and a maximum point a_(r) in the horizontal direction X aswell as a minimum point b_(d) and a maximum point b_(u) in the verticaldirection Y, wherein the four limit points are respectively denoted byO_(l)(a_(l),b_(l)),O_(r) (a_(r),b_(r)),O_(d) (a_(d),b_(d)),O_(u) (a_(u),b_(u)), an inside of a quadrangle formed by the four limit points iscounted as an internal side, and an outside of the quadrangle is countedas an external side; step 3.2. extracting a convex envelope of an axialcenter distribution contour at a minimum slope with the foregoing limitpoints as starting points by anticlockwise traversing all over thepositions of the axial center at all times; and calculating a convexenvelope between a limit point O_(d) and a limit point O_(r) through amethod comprising the following steps; (1) setting a line connecting thepoint O_(d) to the point O_(r) as L₁, wherein a slope α₁ of the line isexpressed as: $\begin{matrix}{\alpha_{1} = \frac{b_{d} - b_{r}}{a_{d} - a_{r}}} & (4)\end{matrix}$ (2) if a set of all common axial center points on anexternal side of the line L₁ is assumed as P={p₁ (a₁,b₁), p₂ (a₂ b₂), p₃(a₃b₃), . . . }, calculating a slope K={β₁,β₂,β₃, . . . } of a lineconnecting the point O_(d) to any point in P; if there are multiplepoints with a same slope, calculating a distance D={dis₁,dis₂,dis₃, . .. } between the corresponding point and the point O_(d); and seeking apoint p′(a_(p),b_(p)) in the set P, and enabling p′ to meet thefollowing formula:β′=minK,β′≤α ₁, and dis′=maxD  (5) wherein, in formula (5), β′represents a slope of a line connecting the point p′ to the point O_(d),and dis′ represents a distance between the point p′ and the point O_(d);and the point p′ is a convex envelope point; (3) replacing the pointO_(d) with the convex envelope point p′, setting a line connecting thepoint p′ to the point O_(r) as L′₁ and a slope of the line L′₁ as α′₁,and then carrying out a next iteration to seek a new convex envelopepoint; (4) repeating step (1), step (2), and step (3); and when adistance between the new convex envelope point and the point O_(r) as 0,stopping the iteration to obtain a convex envelope B′_(dr)={p′₁,p′₂,p′₃,. . . } of an axial center distribution contour between the limit pointO_(d) and the limit point O_(r); and calculating a convex envelopeB′_(ru) between the limit point O_(r) and a limit point O_(u), a convexenvelope B′_(ul) between the limit point O_(u) and a limit point O_(l),and a convex envelope B′_(ld) between the limit point O_(l) and thelimit point O_(d) to obtain a set B_(tu)={O_(d), B′_(dr),O_(r),B′_(ru),O_(u), B′_(ul),O_(l),B′_(ld)} of all convex envelope points in the axialcenter distribution; (5) calculating an area S1 of the convex envelopeof the axial center distribution contour by means of formula (6) below:$\begin{matrix}{{S1} = {\frac{1}{2}{\sum\limits_{e = 1}^{c1}( {{a_{pe}*b_{p{({e + 1})}}} - {a_{p{({e + 1})}}*b_{pe}}} )}}} & (6)\end{matrix}$ wherein, in formula (6), c1 represents the number of theconvex envelope points; step 3.3. calculating a concave envelope of theaxial center distribution according to the convex envelope obtained instep 3.2; and calculating a concave envelope between the limit pointO_(d) and the limit point O_(r) by successively anticlockwise carryingout the following calculation on two continuous convex envelope pointsp′₁(a_(p1),b_(p1)) and p′₂(a_(p2),b_(p2)) in the convex envelopesB′_(dr)={p′₁, p′₂, p′₃, . . . }: (1) setting a line connecting the pointp′₁(a_(p1),b_(p1)) to the point p′₂(a_(p2),b_(p2)) as L₂, wherein aslope α₂ of the line L′ is expressed as: $\begin{matrix}{\alpha_{2} = \frac{b_{p\; 2} - b_{p\; 1}}{a_{p\; 2} - a_{p\; 1}}} & (7)\end{matrix}$ (2) if a set of all common axial center points on aninternal side of the line L₂ is assumed as Q={q₁(a₁, b₁),q₂ (a₂,b₂),q₃(a₃,b₃), . . . }, calculating a slope K′={β′₁,β′₂, β′₃, . . . } ofa line connecting the point p′₁ to any point in Q; if there are multiplepoints with a same slope, calculating a distance D′={dis′₁,dis′₂,dis′₃,. . . } between the corresponding point and the point p′₁; and seeking apoint q′(a_(q),b_(q)) in the set Q, and enabling q′ to meet thefollowing formula:β″=min K′,β″≥α ₂, and dis″=max D′  (8) wherein, in formula (8), β″represents a slope of a line connecting the point q′ to the point p′₁,and dis″ represents a distance between the point q′ and the point p′₁;and the point q′ is a concave envelope point; (3) replacing the pointp′₂ with the concave envelope point q′, setting a line connecting thepoint p′₁ to the point q′ as L′₂ and a slope of the line L′₂ as α₂, andthen carrying out a next iteration to seek a new concave envelope point;(4) repeating step (1), step (2), and step (3); and when a distancebetween the new concave envelope point and the point p′₁ is not greaterthan M, stopping the iteration to obtain a concave envelopeB″_(dr)={q′₁, q′₂, q′₃, . . . } of the axial center distribution contourbetween the limit point O_(d) and the limit point O_(r); wherein, aninitial value of M indicates an average distance between two adjacentaxial center points, and M is calculated by means of formula (9) incombination with formula (10): $\begin{matrix}{\mspace{20mu}{M = \frac{2M^{\prime}}{m( {m - 1} )}}} & (9) \\{M^{\prime} = {{\sum\limits_{i = 2}^{m}{\sqrt{( {a_{1} - a_{i}} )^{2} + ( {b_{1} - b_{i}} )^{2}}}} + {\sum\limits_{i = 3}^{m}{\sqrt{( {a_{2} - a_{i}} )^{2} + ( {b_{1} - b_{i}} )^{2}}}} + \ldots + {\sum\limits_{i = m}^{m}{\sqrt{( {a_{m - 1} - a_{i}} )^{2} + ( {b_{m - 1} - b_{i}} )^{2}}}}}} & (10)\end{matrix}$ calculating a concave envelope B″_(ru) between the limitpoint O_(r) and a limit point O_(u), a concave envelope B″_(ul) betweenthe limit point O_(u) and a limit point O_(l), and a concave envelopeB″_(ld) between the limit point O_(l) and the limit point O_(d) toobtain a setB_(ao)={O_(d),B″_(dr),O_(r),B″_(ru)O_(u),B″_(ul),O_(l),B″_(ld)} of allconcave envelope points in the axial center distribution; (5)calculating an area S2 of the concave envelope of the axial centerdistribution contour by means of formula (11) below: $\begin{matrix}{{S2} = {\frac{1}{2}{\sum\limits_{g = 1}^{c2}( {{a_{qg}*b_{q{({g + 1})}}} - {a_{q{({g + 1})}}*b_{qg}}} )}}} & (11)\end{matrix}$ wherein, in formula (11), c2 represents the number of theconcave envelope points; step 3.4. determining whether or not the setB_(ao), obtained in step 3.3, of the concave envelope points is theenvelope feature of the axial center distribution of the piston rod; (1)calculating a relative error E of S2 and S1 by means of formula (12)till E is less than or equal to 5%, wherein the set B_(ao), obtained instep 3.3, of the concave envelope points is the envelope feature of theaxial center distribution of the piston rod after the calculation isstopped; $\begin{matrix}{{E = {\frac{{{S\; 2} - {S\; 1}}}{S1} \times 100\%}};} & (12)\end{matrix}$ and (2) in a case where the distance M is reduced by 50%when E is greater than 5%, replacing S1 with S2, repeating step 3.3 toobtain a new set B′_(ao) of the concave envelope points as well as anarea S2′ of the concave envelope; and repeatedly calculating a relativeerror E′ of S2′ and S2 by means of formula (13) till E′ is less than orequal to 5%, wherein the iteration is stopped at this moment, and theset B′_(ao) obtained in Step 3.3 is the envelope feature of the axialcenter distribution of the piston rod during the last iteration;$\begin{matrix}{E^{\prime} = {\frac{{{S\; 2^{\prime}} - {S\; 2}}}{S\; 2} \times 100\%}} & (13)\end{matrix}$ step
 4. calculating an information entropy feature of theaxial center distribution O: calculating, by means of formula (14),arithmetic square roots of coordinates of all the points in the axialcenter distribution O to obtain a set S_(m)={s₁ s₂, s₃, s_(m)} of thearithmetic square roots, and then calculating the information entropyfeature Sh of the axial center distribution O by means of formula (15),wherein an initial feature set T={B_(ao),Sh} is formed by theinformation entropy feature Sh and the envelope feature; $\begin{matrix}{{{S_{m}( i^{\prime} )} = {\sqrt{a_{k}^{2} + b_{k}^{2}}}},k,{i^{\prime} = 1},2,\ldots\mspace{14mu},m} & (14) \\{{Sh} = {- {\sum\limits_{i^{\prime} = 1}^{m}( {\frac{{S_{m}( i^{\prime} )}}{\sum\limits_{i^{\prime} = 1}^{m}{{S_{m}( i^{\prime} )}}}*{\log_{2}( \frac{{S_{m}( i^{\prime} )}}{\sum\limits_{i^{\prime} = 1}^{m}{{S_{m}( i^{\prime} )}}} )}} )}}} & (15)\end{matrix}$ step
 5. carrying out, by means of t-distributed stochasticneighbor embedding (T-SNE), an unsupervised dimensionality reduction onthe initial feature set to extract sensitive features of a load; andassuming that the initial feature set T includes 1*Col-dimensionalfeatures, given perplexity is 30, and a given learning rate is 1e-5,setting a label as Labels={0,1,2, . . . ,w},and then inputting theinitial feature set T to a T-SNE algorithm for the unsuperviseddimensionality reduction in the (w+1) load conditions to obtain a setT′={t₁, t₂} of the sensitive features of a 1*2-dimensional load; andstep
 6. firstly, sorting data acquired by the on-line monitoring systemin the (w+1) load conditions to form a training set and a test set;secondly, processing the data in the training set and the test setthrough the above steps to obtain a final training set Train_T′ and afinal test set Test_T′; and thirdly, setting, according to differentreciprocating machinery, the number of neurons of a back-propagation(BP) neural network as 20-30, a learning rate as 0.0005-0.001, trainingaccuracy as 0.0001-0.0005, and maximum iterations as 70-100, theninputting the data set Train_T′ to the BP neural network for training toobtain a classifier capable of distinguishing the (w+1) load conditionsof the reciprocating machinery, and testing the classifier of the BPneural network by means of the test set Test_T′.